function [yc, yinit, objs, times] = cefmc(X, lambda, beta, c, k)

% shapes
V = length(X); 
N = size(X{1}, 1);
times = zeros(1, V+1);


%% Init 
%  --- client ---
% 1. init Hv
H = cell(1, V); dims = zeros(1, V);
for v = 1:V
    tic;
    [uu, ~, ~] = svds(X{v}, c);
    H{v} = uu;
    dims(v) = size(H{v}, 2);
    times(1, v) = toc;
end
% --- server ---
% 2. init yc, Cc, yv, Cv
tic;
[yinit, Cinit] = kmeans(cell2mat(H), k, 'MaxIter',100, 'Replicates',10, 'EmptyAction','drop');
yc = yinit; % Cc = Cinit;
[uu, ~, vv] = svds(Cinit, k);
Cc = uu*vv'.*(k/N); Cc2 = mat2cell(Cc, [k], dims);
times(1, V+1) = toc;
% --- server send to client ----
y = cell(1, V);
for v=1:V
    y{v} = yinit;
%     Y{v} = ind2vec(y{v}', k)'; 
%     [uu, ~, vv] = svds(H{v}'*Y{v}, k);
%     C{v} = vv*uu'.*(c/N);
end
C = mat2cell(Cinit, [k], dims);


%% Clustering
max_iter = 100; in_iter = 5;
Isconverg = 0; t = 1; obj = 0;
while(Isconverg == 0)
    
    % --- client ---
    
    for v = 1:V
        tic;
        % 1. update Hv
        for i = 1:in_iter
            M{v} = X{v}*(X{v}'*H{v}) + lambda*C{v}(y{v},:);
            [uu, ~, vv] = svds(M{v}, c);
            H{v} = uu*vv';
        end

%         % 2. update Cv
%         Y{v} = ind2vec(y{v}', k)'; 
%         [uu, ~, vv] = svds(H{v}'*Y{v}+Cc2{v}(yc,:)'*Y{v}, k);
%         C{v} = vv*uu'.*(c/N);

         % 3. update Yv
        S{v} = lambda*H{v}*C{v}' + beta*Cc2{v}(yc,:)*C{v}';
        [~, tmp] = max(S{v}, [], 2);
        y{v} = tmp;
        
        times(t+1, v) = toc;
    end
    
    % --- server ---
    tic;
    YCcon = cell(1, V);
    for v = 1:V
        YCcon{v} = C{v}(y{v},:);
    end
    YCcon = cell2mat(YCcon);
    
    % 4. update Cc
    Yc = ind2vec(yc', k)'; 
    [uu, ~, vv] = svds(YCcon'*Yc, k);
    Cc = vv*uu'.*(k/N); Cc2 = mat2cell(Cc, [k], dims);
    
    % 5. update Yc
    [~, tmp] = max(YCcon*Cc', [], 2);
    yc = tmp; 
    
    times(t+1, V+1) = toc;
    
    % compute objective value
    obj1 = 0; obj2 = 0; obj3 = 0;
    for v = 1:V
        tic;
        obj1 = obj1 + trace((H{v}'*X{v})*(X{v}'*H{v}));
        obj2 = obj2 + trace(H{v}'*C{v}(y{v},:));
        times(t+1, v) = times(t+1, v) + toc;
        tic;
        obj3 = obj3 + trace(C{v}(y{v},:)'*Cc2{v}(yc,:));
        times(t+1, V+1) = times(t+1, V+1) + toc;
    end
    obj(t) = obj1 + lambda*obj2 + beta*obj3;
    objs{1}(t) = obj(t); objs{2}(t) = obj1; objs{3}(t) = obj2; objs{4}(t) = obj3;
    
    % if convergent
    if t>1 && (abs((obj(t-1)-obj(t))/obj(t))<1e-5 || t>max_iter)
        Isconverg = 1;
    end
    
    t = t + 1;
    
end


end

